Heritability and genetic correlates of childhood behavioural traits

excl <- c(
  "IcqFull","CbclFull","EasFull","CfqFull","ItseaFull","ScqFull","SdqFull",
  "CccFull","SprakFull","NhpicFull","RsdbdFull","PeqFull","SbFull",
  "SclFull","SppaFull","CbclAdhd"
)

SsSetM <- "ChPsyGenFscAsr"

h2TabM <- 
  fread(paste0("Tables_GenSem/ChPsy/mvLDSC_",SsSetM,"_h2FullTable.txt"),header=T) %>% 
  # select(Trait=V1,h2=V2,SE=V3,Int=V4,IntSE=V5,Ratio=V6,RatioSE=V7,Z=V8) %>% 
  # select(Trait, h2, SE, Z)  %>% 
  # mutate(across(c(h2, SE, Z), ~round(.x,2)))  %>% 
  distinct(Trait,.keep_all = T) %>% 
  mutate(Trait=str_replace_all(Trait,"Q0","_Q0")) %>% 
  mutate(Trait=str_replace_all(Trait,"Q1","_Q1")) %>% 
  separate(Trait,c("Trait","Age"),"_") %>% 
  mutate(alpha = ifelse(h2<0 | h2-SE<0, 1, 1)) %>% 
  mutate(Type = paste("GWA")) %>% 
  mutate(Trait=str_replace_all(Trait,"Ch","")) %>% 
  mutate(Age=str_replace_all(Age,"Q04","6mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q05","18mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q06","3yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q07","5yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q09","8yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q11C","14yr(c)")) %>% 
  mutate(Age=str_replace_all(Age,"Q11M","14yr(m)")) %>% 
  filter(!Trait %in% excl)


# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SsSetV <- "ChPsyGenQrsAsr"

h2TabV <- 
  fread(paste0("Tables_GenSem/ChPsy/mvLDSC_",SsSetV,"_h2FullTable.txt"),header=T) %>% #,sep2=" ") %>% 
  # select(Trait=V1,h2=V2,SE=V3,Int=V4,IntSE=V5,Ratio=V6,RatioSE=V7,Z=V8) %>% 
  # select(Trait, h2, SE, Z)  %>% 
  # mutate(across(c(h2, SE, Z), ~round(.x,2)))  %>% 
  distinct(Trait,.keep_all = T) %>% 
  mutate(Trait=str_replace_all(Trait,"Q0","_Q0")) %>% 
  mutate(Trait=str_replace_all(Trait,"Q1","_Q1")) %>% 
  separate(Trait,c("Trait","Age")) %>% 
  mutate(alpha = ifelse(h2<0 | h2-SE<0, 1, 1))%>% 
  mutate(Type = paste("vGWA")) %>% 
  mutate(Trait=str_replace_all(Trait,"Ch","")) %>% 
  mutate(Age=str_replace_all(Age,"Q04","6mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q05","18mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q06","3yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q07","5yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q09","8yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q11C","14yr(c)")) %>% 
  mutate(Age=str_replace_all(Age,"Q11M","14yr(m)")) %>% 
  filter(!Trait %in% excl)

h2Tab1 <- rbind(h2TabM,h2TabV)

h2Tab2 <- 
  h2TabM %>% 
  select(Trait,Age,h2_M=h2,se_M=SE,Z_M=Z) %>% 
  full_join(h2TabV %>% select(Trait,Age,h2_V=h2,se_V=SE,Z_V=Z))


N <- fread("Tables_GenSem/ChPsy/ChPsyGenQrsAsr_SampleSizes.txt",header=F) %>% 
  mutate(V1=str_replace_all(V1,"_RegenieMQt.Munged.sumstats.gz","")) %>% 
  mutate(V1=str_replace_all(V1,"_RegenieVQt.Munged.sumstats.gz","")) %>% 
  mutate(V1=str_replace_all(V1,"GwaSumStats/ChPsyGenQrsAsr/","")) %>% 
  # mutate(V1=str_replace_all(V1,"([a-z])([A-Z])", "\\1 \\2")) %>% 
  mutate(V1=str_replace_all(V1,"Q", "_Q")) %>% 
  separate(V1, c("Trait","Age")) %>% 
  select(Trait, Age, N=V2) %>% 
  mutate(Trait=str_replace_all(Trait,"Ch","")) %>% 
  mutate(Age=str_replace_all(Age,"Q04","6mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q05","18mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q06","3yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q07","5yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q09","8yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q11C","14yr(c)")) %>% 
  mutate(Age=str_replace_all(Age,"Q11M","14yr(m)")) %>% 
  filter(!Trait %in% excl)

ModStat <- fread("Tables_Phe/xIRT_ModelStats.tsv") %>%
  mutate(Trait=paste0(Scale,Subscale)) %>% 
  select(Trait,Age=Wave,everything()) %>% 
  mutate(Trait=str_replace_all(Trait,"Ch","")) %>% 
  mutate(Age=str_replace_all(Age,"Q04","6mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q05","18mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q06","3yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q07","5yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q09","8yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q11C","14yr(c)")) %>% 
  mutate(Age=str_replace_all(Age,"Q11M","14yr(m)")) %>% 
  select(!matches("Scale"),!matches("Subscale")) %>% 
  inner_join(h2Tab2) %>% inner_join(N) %>% 
  filter(!Trait %in% excl)

ModStat$Age <- factor(ModStat$Age,
                         levels=c("6mo","18mo","3yr","5yr","8yr","14yr(m)","14yr(c)"),
                         labels=c("6 m/o","18 m/o","3 y/o","5 y/o","8 y/o","14 y/o","14 y/o")
                         )


# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

mCorMat <- fread("Tables_GenSem/ChPsy/MvLDSC_ChPsyGenFscAsr_rG_est.csv") %>% as.matrix()
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Ch","")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Q04","_6mo")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Q05","_18mo")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Q06","_3yr")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Q07","_5yr")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Q09","_8yr")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat), "Q11C","_14yr(c)")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat), "Q11M","_14yr(m)")
row.names(mCorMat) <- colnames(mCorMat)


mCorTab <- round(mCorMat,2)
# mCorTab[upper.tri(mCorTab)] <- NA
mCorTab <- reshape2::melt(mCorTab)
names(mCorTab)[3] <- "rg"

ErrMat  <- fread("Tables_GenSem/ChPsy/MvLDSC_ChPsyGenFscAsr_rG_err.csv") %>% as.matrix()
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Ch","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q04","_6mo")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q05","_18mo")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q06","_3yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q07","_5yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q09","_8yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q11C","_14yr(c)")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q11M","_14yr(m)")
row.names(ErrMat) <- colnames(ErrMat)

ErrTab  <- round(ErrMat,3)
# ErrTab[upper.tri(ErrTab)] <- NA
ErrTab  <- reshape2::melt(ErrTab)
names(ErrTab)[3] <- "SE"

mCorTab <- mCorTab %>% full_join(ErrTab,by=c("Var1","Var2"))

mCorTab$Z <- mCorTab$rg/mCorTab$SE

Pval <- fread("Tables_GenSem/ChPsy/mvLDSC_ChPsyGenFscAsr_rG_pval.txt")%>% select(Var1=V1,Var2=V2,Pval=V3) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"_RegenieVQt.Munged.sumstats.gz",""))%>% 
  mutate(across(c(Var1,Var2),str_replace_all,"_RegenieMQt.Munged.sumstats.gz",""))%>% 
  mutate(across(c(Var1,Var2),str_remove_all,"Ch")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q04","_6mo")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q05","_18mo")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q06","_3yr")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q07","_5yr")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q09","_8yr")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q11C","_14yr(c)")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q11M","_14yr(m)")) %>% 
  mutate(across(c(Var1,Var2),as.factor))
  
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Ch","")

mCorTab <- mCorTab %>% full_join(Pval,by=c("Var1","Var2")) %>% na.omit()  %>% 
  filter(!grepl(x=Var1,pattern=paste0(excl,collapse="|"))) %>% 
  filter(!grepl(x=Var2,pattern=paste0(excl,collapse="|"))) %>% 
  mutate(Z=round(Z,1))


# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

vCorMat <- fread("Tables_GenSem/ChPsy/MvLDSC_ChPsyGenQrsAsr_rG_est.csv") %>% as.matrix()
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"_RegenieVQt.Munged.sumstats.gz","")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Ch","")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q04","6mo")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q05","18mo")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q06","3yr")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q07","5yr")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q09","8yr")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q11C","14yr(c)")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q11M","14yr(m)")
row.names(vCorMat) <- colnames(vCorMat)

vCorTab <- round(vCorMat,2)
# vCorTab[upper.tri(vCorTab)] <- NA
vCorTab <- reshape2::melt(vCorTab)
names(vCorTab)[3] <- "rg"



ErrMat  <- fread("Tables_GenSem/ChPsy/MvLDSC_ChPsyGenQrsAsr_rG_err.csv") %>% as.matrix()
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"_RegenieVQt.Munged.sumstats.gz","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Ch","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q04","6mo")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q05","18mo")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q06","3yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q07","5yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q09","8yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q11C","14yr(c)")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q11M","14yr(m)")
row.names(ErrMat) <- colnames(ErrMat)

ErrTab  <- round(ErrMat,3)
# ErrTab[upper.tri(ErrTab)] <- NA
ErrTab  <- reshape2::melt(ErrTab)
names(ErrTab)[3] <- "SE"

vCorTab <- vCorTab %>% full_join(ErrTab,by=c("Var1","Var2"))

vCorTab$Z <- vCorTab$rg/vCorTab$SE

Pval <- fread("Tables_GenSem/ChPsy/mvLDSC_ChPsyGenQrsAsr_rG_pval.txt")%>% 
  select(Var1=V1,Var2=V2,Pval=V3) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"_RegenieVQt.Munged.sumstats.gz",""))%>% 
  mutate(across(c(Var1,Var2),str_replace_all,"_RegenieMQt.Munged.sumstats.gz",""))%>% 
  mutate(across(c(Var1,Var2),as.factor)) %>% 
  mutate(across(c(Var1,Var2),str_remove_all,"Ch")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q04","6mo")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q05","18mo")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q06","3yr")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q07","5yr")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q09","8yr")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q11C","14yr(c)")) %>% 
  mutate(across(c(Var1,Var2),str_replace_all,"Q11M","14yr(m)")) %>% 
  mutate(across(c(Var1,Var2),as.factor))

vCorTab <- vCorTab %>% full_join(Pval,by=c("Var1","Var2")) %>% na.omit()%>% 
  mutate(across(c(Var1,Var2),as.character))%>% 
  filter(!grepl(x=Var1,pattern=paste0(excl,collapse="|"))) %>% 
  filter(!grepl(x=Var2,pattern=paste0(excl,collapse="|"))) %>% 
  mutate(Z=round(Z,1))


meanh2V  <- round(mean(h2TabV$h2),2)
meanh2m  <- round(mean(h2TabM$h2),2)
meanZV   <- round(mean(h2TabV$Z),2)
meanZm   <- round(mean(h2TabM$Z),2)
meanh2Vz <- round(mean(h2TabV %>% filter(Z>3) %>% .$h2),2)
meanh2mz <- round(mean(h2TabM %>% filter(Z>3) %>% .$h2),2)
h2V_N  <- nrow(h2TabV %>% filter(Z > 3))
h2M_N  <- nrow(h2TabM %>% filter(Z > 3))
meanXV   <- round(mean(h2TabV$MeanChi2),2)
meanXm   <- round(mean(h2TabM$MeanChi2),2)
# max(h2TabV$h2)
# max(h2TabM$h2)
RgV_N  <- nrow(vCorTab %>% filter(Z > 3))
RgM_N  <- nrow(mCorTab %>% filter(Z > 3))

meanVRgZ  <- round(mean(abs(vCorTab$Z)),2)
meanMRgZ  <- round(mean(abs(mCorTab$Z)),2)

Heritability

Tables

h2Tab2$Age <- factor(h2Tab2$Age,levels=c("6mo","18mo","3yr","5yr","8yr","14yr(m)","14yr(c)"))

datatable(
  h2Tab2 %>% purrr::set_names(c("Trait","Age","h2 [GWA]","SE [GWA]","Z [GWA]","h2 [vGWA]","SE [vGWA]","Z [vGWA]")), # %>% mutate(`|Zm|` = abs()),
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

datatable(
  h2Tab1, # %>% mutate(`|Zm|` = abs()),
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

Mean h2-vGWA = 0.02, meanZ= 1.22; Mean h2-mGWA = 0.06, meanZ= 3.76;
when Z>3: Mean h2-vGWA = 0.05; Mean h2-mGWA = 0.08

Plots

h2 (I)

ModStat <- fread("Tables_Phe/xIRT_ModelStats.tsv") %>%
  mutate(Trait=paste0(Scale,Subscale)) %>% 
  select(Trait,Age=Wave,everything()) %>% 
  mutate(Trait=str_replace_all(Trait,"Ch","")) %>% 
  mutate(Age=str_replace_all(Age,"Q04","6mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q05","18mo")) %>% 
  mutate(Age=str_replace_all(Age,"Q06","3yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q07","5yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q09","8yr")) %>% 
  mutate(Age=str_replace_all(Age,"Q11C","14yr(c)")) %>% 
  mutate(Age=str_replace_all(Age,"Q11M","14yr(m)")) %>% 
  select(!matches("Scale"),!matches("Subscale")) %>% 
  inner_join(h2Tab2) %>% inner_join(N) %>% 
  filter(!Trait %in% excl) %>% 
  mutate(Age=factor(Age,levels=c("6mo","18mo","3yr","5yr","8yr","14yr(m)","14yr(c)"),
                        labels=c("6 m/o","18 m/o","3 y/o","5 y/o","8 y/o","14 y/o","14 y/o")
                         )) %>% 
  mutate(
    h2_GWA=round(h2_M,2),h2_vGWA=round(h2_V,2),
    se_GWA=round(se_M,2),se_vGWA=round(se_V,2),
    lse_vGWA=round(h2_V-se_V,2),use_vGWA=round(h2_V+se_V,2),
    lse_GWA=round(h2_M-se_M,2), use_GWA=round(h2_M+se_M,2)
         )
    
Mh2_Vh2 <-
  ggplot(
    ModStat,
    aes(
      # label=Trait,label2=h2_GWA,label3=se_GWA,label4=h2_vGWA,label5=se_vGWA #, #))+
      x=h2_GWA,y=h2_vGWA
      )) +
  geom_linerange(aes(ymin=round(h2_V-se_V,2),ymax=round(h2_V+se_V,2),color=Trait)) +
  geom_linerange(aes(xmin=round(h2_M-se_M,2),xmax=round(h2_M+se_M,2),color=Trait)) +
  # scale_x_continuous(limits=c(-0.06,0.16),breaks=c(seq(0,0.15,0.05)))+
  # scale_y_continuous(limits=c(-0.06,0.16),breaks=c(seq(0,0.15,0.05)))+
  scale_x_continuous(limits=c(-0.1,0.2),breaks=c(seq(-.1,.2,0.1)))+
  scale_y_continuous(limits=c(-0.1,0.2),breaks=c(seq(-.1,.2,0.1)))+
  geom_hline(yintercept=0,colour="limegreen",linetype="dotted") +
  geom_vline(xintercept=0,colour="limegreen",linetype="dotted") +
  # geom_smooth(method=lm,se=FALSE,colour="limegreen") +
  facet_wrap(facets=vars(Age),ncol=3) + #, scales = "free",space = "free")+
  theme_minimal() +
  scale_colour_viridis_d(option="plasma") +
  guides(color="none") #+
  # xlab(bquote({h^2}[LDSC]))+ 
  # ylab(bquote({h^2}[vLDSC]))

Mh2_Vh2_ly <- plotly::ggplotly(Mh2_Vh2, tooltip = c("Trait","h2_GWA","se_GWA","h2_vGWA","se_vGWA"))
saveWidget(Mh2_Vh2_ly, "./ChPsy_h2Plotly.html", selfcontained = T, libdir = "lib")

Mh2_Vh2_ly

h2 (II)

h2Tab1$Age <- factor(h2Tab1$Age,levels=c("6mo","18mo","3yr","5yr","8yr","14yr(m)","14yr(c)"))

scaleFUN <- function(x) sprintf("%.2f", x)

h2Plot<-
  ggplot(h2Tab1 %>% mutate(h2=round(h2,2)), 
       aes(x=Type, y = h2, colour = Type, alpha = alpha )) +
  geom_point(aes(size = Z)) +
  geom_errorbar(aes(ymin = h2 - SE, ymax = h2 + SE), width = 0.2) +
  facet_grid(Trait~Age,scales = "free_y",) +
  theme_minimal()+
  theme(
    legend.position="top", 
    legend.title = element_blank(),
    axis.text.x=element_blank(),
    axis.title.x=element_blank(),
    strip.text.y=element_text(angle=0)
    )+ 
  scale_colour_viridis_d(option="plasma",begin=0.5,end=0.8)+
  guides(size="none",alpha="none")+
  geom_hline(yintercept=0,colour="limegreen",linetype="dashed")+
  scale_y_continuous(labels=scaleFUN)

h2Plot

Reliability

Rel_h2 <-
  ggplot(
    ModStat %>% 
      select(Trait, vGWA=h2_V, mGWA=h2_M, Reliability) %>% 
      pivot_longer(!c(Trait,Reliability),names_to="Effect",values_to="h2"),
    aes(x=Reliability, y=h2, group=Effect, colour=Effect,shape=Effect,label=Trait)) +
  geom_point() +
  geom_smooth(method=lm,se=FALSE) +
  scale_colour_viridis_d(option="viridis",direction=-1,begin=0.4,end=0.8)+
  theme_minimal() +
  directlabels::geom_dl(aes(label=Effect),method="smart.grid")+
  guides(color="none",shape="none")#+ 
  # ylab(bquote({h^2}))


Rel_h2_ly <- plotly::ggplotly(Rel_h2, tooltip = c("Reliability","h2","Effect","Trait"))
saveWidget(Mh2_Vh2_ly, "./ChPsy_h2Plotly.html", selfcontained = T, libdir = "lib")

Rel_h2_ly

# ModCor<-polycor::hetcor(ModStat %>% select(LL,AIC,BIC,M2,df,RMSEA,SRMSR,TLI,CFI,Reliability,h2_M,Z_M,h2_V,Z_V,N),ML=T,use="p")

# datatable(cortably(ModCor),
#   filter='bottom',extensions='Buttons',rownames=F,
#   options=list(
#     dom='lfrtiBp',buttons=c('csv'),
#     searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
#     columnDefs=list(list(className='dt-left',targets="_all"))
#     ))

Genetic correlations

Tables

Factor scores (GWA)

datatable(
  mCorTab %>% select(Var1,Var2,rg,SE,Pval,Z),
  filter='bottom',rownames=F,extensions='Buttons',
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

# max(abs(mCorTab$Z))
# mean(abs(mCorTab$Z))
# max(abs(vCorTab$Z))
# mean(abs(vCorTab$Z))

Quantile rank scores (vGWA)

datatable(
  vCorTab %>% select(Var1,Var2,rg,SE,Pval,Z),
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

Plots

Factor scores (GWA)

Full

mCor  <- mCorTab %>% select(Var1,Var2,rg,SE,Z)
vars  <- unique(c(levels(mCor$Var1),levels(mCor$Var2)))

excl <- c(
  "IcqFull","CbclFull","EasFull","CfqFull","ItseaFull","ScqFull","SdqFull",
  "CccFull","SprakFull","NhpicFull","RsdbdFull","PeqFull","SbFull",
  "SclFull","SppaFull","CbclAdhd"
)

vars  <- vars[!grepl(paste0(excl,collapse="|"),vars)]

combos <- apply(combn(vars,2),1,paste) %>% as.data.frame()
autos  <- data.frame(V1=vars,V2=vars) 
combos <- combos %>% bind_rows(autos) %>% arrange(V1,V2)
names(combos)<-c("Var1","Var2")

mCor      <- 
  bind_rows(
    mCorTab %>% select(Var1,Var2,rg),
    mCorTab %>% select(Var1=Var2,Var2=Var1,rg)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(rg=ifelse(Var1==Var2,1,rg)) %>% 
  pivot_wider(names_from=Var1,values_from=rg,values_fill=NA) %>%
  column_to_rownames("Var2") %>% 
  as.matrix()

mCor[is.na(mCor)] <- 0

mCor[mCor>1] <- 1
mCor[mCor<(-1)] <- (-1)

Err  <- 
  bind_rows(
    mCorTab %>% select(Var1,Var2,SE),
    mCorTab %>% select(Var1=Var2,Var2=Var1,SE)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(SE=ifelse(Var1==Var2,0,SE)) %>% 
 pivot_wider(names_from=Var1,values_from=SE) %>%
 column_to_rownames("Var2") %>% 
 as.matrix()

Err[is.na(Err)] <- 0

Zsc  <- 
  bind_rows(
    mCorTab %>% select(Var1,Var2,Z),
    mCorTab %>% select(Var1=Var2,Var2=Var1,Z)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(Z=ifelse(Var1==Var2,2,Z)) %>% 
 pivot_wider(names_from=Var1,values_from=Z) %>%
 column_to_rownames("Var2") %>% 
 as.matrix()

Zsc[is.na(Zsc)] <- 0
Zsc  <- abs(Zsc)


ChmCorplotly<-
    heatmaply::heatmaply_cor(
      mCor,
      node_type = "scatter",
      # point_size_mat = -log(Pval),
      # point_size_name = "-logP",
      point_size_mat=Zsc,
      point_size_name="|Z|",
      label_names = c("x", "y", "rg"),
      fontsize_col=8,
      fontsize_row=8,
      grid_size = 0.001,
      dendrogram = "both",
      subplot_heights=c(0.06, 0.94),
      subplot_widths=c(0.94,0.06),
      hide_colorbar=T
    ) 

ChmCorplotly$height <- nrow(mCor)*15
ChmCorplotly$width  <- nrow(mCor)*15
ChmCorplotly$x$data[[4]]$marker$size <- (ChmCorplotly$x$data[[4]]$marker$size/2)

ChmCorplotly

|Z| > 3

mCor   <- mCorTab %>% select(Var1,Var2,rg,SE,Z)
vars   <- unique(c(levels(mCor$Var1),levels(mCor$Var2)))
vars   <- vars[!grepl(paste0(excl,collapse="|"),vars)]
combos <- apply(combn(vars,2),1,paste) %>% as.data.frame()
autos  <- data.frame(V1=vars,V2=vars) 
combos <- combos %>% bind_rows(autos) %>% arrange(V1,V2)
names(combos)<-c("Var1","Var2")

mCor      <- 
  bind_rows(
    mCorTab %>% filter(abs(Z)>3) %>%select(Var1,Var2,rg),
    mCorTab %>% filter(abs(Z)>3) %>%select(Var1=Var2,Var2=Var1,rg)
  ) %>%
  inner_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(rg=ifelse(Var1==Var2,1,rg)) %>% 
  pivot_wider(names_from=Var1,values_from=rg,values_fill=NA) %>%
  column_to_rownames("Var2") %>% 
  as.matrix()

mCor[is.na(mCor)] <- 0

mCor[mCor>1] <- 1
mCor[mCor<(-1)] <- (-1)

Err  <- 
  bind_rows(
    mCorTab %>% filter(abs(Z)>3) %>%select(Var1,Var2,SE),
    mCorTab %>% filter(abs(Z)>3) %>%select(Var1=Var2,Var2=Var1,SE)
  ) %>%
  inner_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(SE=ifelse(Var1==Var2,0,SE)) %>% 
 pivot_wider(names_from=Var1,values_from=SE) %>%
 column_to_rownames("Var2") %>% 
 as.matrix()

Err[is.na(Err)] <- 0

Zsc  <- 
  bind_rows(
    mCorTab %>% filter(abs(Z)>3) %>%select(Var1,Var2,Z),
    mCorTab %>% filter(abs(Z)>3) %>%select(Var1=Var2,Var2=Var1,Z)
  ) %>%
  inner_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(Z=ifelse(Var1==Var2,2,Z)) %>% 
 pivot_wider(names_from=Var1,values_from=Z) %>%
 column_to_rownames("Var2") %>% 
 as.matrix()

Zsc[is.na(Zsc)] <- 0
Zsc  <- abs(Zsc)


ChmCorplotly<-
    heatmaply::heatmaply_cor(
      as.data.frame(mCor),
      node_type = "scatter",
      # point_size_mat = -log(Pval),
      # point_size_name = "-logP",
      point_size_mat=Zsc,
      point_size_name="|Z|",
      label_names = c("x", "y", "rg"),
      fontsize_col=8,
      fontsize_row=8,
      grid_size = 0.001,
      dendrogram = "both",
      subplot_heights=c(0.06, 0.94),
      subplot_widths=c(0.94,0.06),
      hide_colorbar=T
    ) 

ChmCorplotly$height <- nrow(mCor)*20
ChmCorplotly$width  <- nrow(mCor)*20
ChmCorplotly$x$data[[4]]$marker$size <- (ChmCorplotly$x$data[[4]]$marker$size/2)

ChmCorplotly

Quantile rank scores (vGWA)

Full

vCor  <- vCorTab %>% select(Var1,Var2,rg,SE,Z)
vars  <- unique(c(vCor$Var1,vCor$Var2))

excl <- c(
  "IcqFull","CbclFull","EasFull","CfqFull","ItseaFull","ScqFull","SdqFull",
  "CccFull","SprakFull","NhpicFull","RsdbdFull","PeqFull","SbFull",
  "SclFull","SppaFull","CbclAdhd"
)

vars  <- vars[!grepl(paste0(excl,collapse="|"),vars)]

combos <- apply(combn(vars,2),1,paste) %>% as.data.frame()
autos  <- data.frame(V1=vars,V2=vars) 
combos <- combos %>% bind_rows(autos) %>% arrange(V1,V2)
names(combos)<-c("Var1","Var2")

vCor      <- 
  bind_rows(
    vCorTab %>% select(Var1,Var2,rg),
    vCorTab %>% select(Var1=Var2,Var2=Var1,rg)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(rg=ifelse(Var1==Var2,1,rg)) %>% 
  pivot_wider(names_from=Var1,values_from=rg,values_fill=NA) %>%
  column_to_rownames("Var2") %>% 
  as.matrix()

vCor[is.na(vCor)] <- 0

vCor[vCor>1] <- 1
vCor[vCor<(-1)] <- (-1)

Err  <- 
  bind_rows(
    vCorTab %>% select(Var1,Var2,SE),
    vCorTab %>% select(Var1=Var2,Var2=Var1,SE)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(SE=ifelse(Var1==Var2,0,SE)) %>% 
 pivot_wider(names_from=Var1,values_from=SE) %>%
 column_to_rownames("Var2") %>% 
 as.matrix()

Err[is.na(Err)] <- 0

Zsc  <- 
  bind_rows(
    vCorTab %>% select(Var1,Var2,Z),
    vCorTab %>% select(Var1=Var2,Var2=Var1,Z)
  ) %>%
  full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(Z=ifelse(Var1==Var2,2,Z)) %>% 
 pivot_wider(names_from=Var1,values_from=Z) %>%
 column_to_rownames("Var2") %>% 
 as.matrix()

Zsc[is.na(Zsc)] <- 0
Zsc  <- abs(Zsc)

ChvCorplotly<-
    heatmaply::heatmaply_cor(
      as.data.frame(vCor),
      node_type = "scatter",
      # point_size_mat = -log(Pval),
      # point_size_name = "-logP",
      point_size_mat=Zsc,
      point_size_name="|Z|",
      label_names = c("x", "y", "rg"),
      fontsize_col=8,
      fontsize_row=8,
      grid_size = 0.001,
      dendrogram = "both",
      subplot_heights=c(0.06, 0.94),
      subplot_widths=c(0.94,0.06),
      hide_colorbar=T
    ) 

ChvCorplotly$height <- nrow(vCor)*20
ChvCorplotly$width  <- nrow(vCor)*20
ChvCorplotly$x$data[[4]]$marker$size <- (ChvCorplotly$x$data[[4]]$marker$size/2)

ChvCorplotly

|Z| > 3

vCor      <- 
  bind_rows(
    vCorTab %>% filter(abs(Z)>3) %>%select(Var1,Var2,rg),
    vCorTab %>% filter(abs(Z)>3) %>%select(Var1=Var2,Var2=Var1,rg)
  ) %>%
  inner_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(rg=ifelse(Var1==Var2,1,rg)) %>% 
  pivot_wider(names_from=Var1,values_from=rg,values_fill=NA) %>%
  column_to_rownames("Var2") %>% 
  as.matrix()

vCor[is.na(vCor)] <- 0

vCor[vCor>1] <- 1
vCor[vCor<(-1)] <- (-1)

Err  <- 
  bind_rows(
    vCorTab %>% filter(abs(Z)>3) %>%select(Var1,Var2,SE),
    vCorTab %>% filter(abs(Z)>3) %>%select(Var1=Var2,Var2=Var1,SE)
  ) %>%
  inner_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(SE=ifelse(Var1==Var2,0,SE)) %>% 
 pivot_wider(names_from=Var1,values_from=SE) %>%
 column_to_rownames("Var2") %>% 
 as.matrix()

Err[is.na(Err)] <- 0

Zsc  <- 
  bind_rows(
    vCorTab %>% filter(abs(Z)>3) %>%select(Var1,Var2,Z),
    vCorTab %>% filter(abs(Z)>3) %>%select(Var1=Var2,Var2=Var1,Z)
  ) %>%
  inner_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>% 
  mutate(Z=ifelse(Var1==Var2,2,Z)) %>% 
 pivot_wider(names_from=Var1,values_from=Z) %>%
 column_to_rownames("Var2") %>% 
 as.matrix()

Zsc[is.na(Zsc)] <- 0
Zsc  <- abs(Zsc)


ChvCorplotly<-
    heatmaply::heatmaply_cor(
      as.data.frame(abs(vCor)),
      node_type = "scatter",
      # point_size_mat = -log(Pval),
      # point_size_name = "-logP",
      point_size_mat=Zsc,
      point_size_name="|Z|",
      label_names = c("x", "y", "rg"),
      fontsize_col=8,
      fontsize_row=8,
      grid_size = 0.001,
      dendrogram = "both",
      subplot_heights=c(0.06, 0.94),
      subplot_widths=c(0.94,0.06),
      hide_colorbar=T
    ) 

ChmCorplotly$height <- nrow(mCor)*20
ChmCorplotly$width  <- nrow(mCor)*20
ChvCorplotly$x$data[[4]]$marker$size <- (ChvCorplotly$x$data[[4]]$marker$size/1.5)

ChvCorplotly

Top SNPs

SnpTab <- fread("Tables_Snps/ChPsy/ClumpedTopSnps.tsv")

SnpTab$af <- round(SnpTab$af,2)
SnpTab$mP <- formatC(SnpTab$mP, format = "e", digits = 2)
SnpTab$vP <- formatC(SnpTab$vP, format = "e", digits = 2)

DT::datatable(
  SnpTab,
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX=T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

cor.test(as.numeric(SnpTab$mP),as.numeric(SnpTab$vP))

## 
##  Pearson's product-moment correlation
## 
## data:  as.numeric(SnpTab$mP) and as.numeric(SnpTab$vP)
## t = -0.36345, df = 384, p-value = 0.7165
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.11814149  0.08142252
## sample estimates:
##         cor 
## -0.01854418

PheWas

Table

vPheWasTab <- fread("Tables_Snps/ChPsy/vPheWasTab.tsv")
mPheWasTab <- fread("Tables_Snps/ChPsy/mPheWasTab.tsv")
PheWasTab3 <- vPheWasTab %>% bind_rows(mPheWasTab)

PheWasLab <- 
  PheWasTab3 %>% 
  arrange(PheWasP) %>% distinct(PheWasTrait,MoBaTrait,.keep_all=T) %>% 
  select(PheWasTrait,MoBaTrait,PheWasP) %>%
  mutate(Label=str_replace_all(PheWasTrait, "\\s*\\([^\\)]+\\)", "")) %>% 
  mutate(Label=str_replace_all(Label, "[[:punct:]]", "")) %>% 
  mutate(Label=str_to_title(Label)) %>% 
  mutate(Label=str_replace_all(Label, " ", "")) %>% 
  mutate(Label=str_replace_all(Label, "CellCount", "Count")) %>% 
  distinct(Label,MoBaTrait,.keep_all=T) %>% 
  inner_join(PheWasTab3,c("MoBaTrait","PheWasTrait","PheWasP")) %>% 
  select(Label,PheWasTrait,MoBaTrait,vP,mP,PheWasP,rsid) %>% distinct() %>% 
  mutate(MoBaTrait=str_replace_all(MoBaTrait, "C_", " ")) %>% 
  mutate(MoBaTrait=str_replace_all(MoBaTrait, "Full", "")) %>% 
  # mutate(MoBaTrait=str_replace_all(MoBaTrait, "([a-z])([A-Z])", "\\1 \\2")) %>% 
  mutate(MoBaTrait=str_replace_all(MoBaTrait, "([a-z])([0-9])", "\\1 \\2"))%>% 
  separate(MoBaTrait,c("MoBaTrait","Age")," ")%>% 
  mutate(Age = factor(Age, levels=c("6mo","18mo","3yr","5yr","8yr","14yr"))) #%>% 
  # mutate(MoBaTrait=paste0(Scale,Trait))

PheWasLab <- 
  bind_rows(
  PheWasLab %>% filter(vP<5e-8) %>% rename("GwaP"=vP) %>% mutate(Effect="vQTL"),
  PheWasLab %>% filter(mP<5e-8) %>% rename("GwaP"=mP) %>% mutate(Effect="mQTL")
  ) %>% 
  mutate(GwaP=round(-log10(GwaP),1),PheWasP=round(-log10(PheWasP),1))


PheWasTab <- PheWasTab3 %>% select(rsid,PheWasTrait,MoBaTrait,PheWasP,mP,vP)

PheWasTab$PheWasP <- formatC(PheWasTab$PheWasP, format = "e", digits = 2)
PheWasTab$mP <- formatC(PheWasTab$mP, format = "e", digits = 2)
PheWasTab$vP <- formatC(PheWasTab$vP, format = "e", digits = 2)

DT::datatable(
  PheWasTab,
  filter='bottom',extensions='Buttons',rownames=F,
  options=list(
    dom='lfrtiBp',buttons=c('csv'),
    searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX=T,
    columnDefs=list(list(className='dt-left',targets="_all"))
    ))

Plots

library(ggplot2)
library(ggrepel)
library(directlabels)
library(dplyr)
library(stringr)
library(data.table)
library(directlabels)

## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  

PheWasPlot <- 
  ggplot(
    PheWasLab,aes(y=PheWasP,x=GwaP,shape=MoBaTrait,colour=MoBaTrait,label=PheWasTrait))+ 
    geom_point()+theme_minimal()+
    facet_grid(cols=vars(Age),rows=vars(Effect))+
    scale_colour_viridis_d(begin=.2,end=.8)+
    scale_shape_manual(values=c(0:14))+
    theme(legend.position="bottom")
  
plotly::ggplotly(PheWasPlot)

# PheWasPlotly <- plotly::ggplotly(PheWasPlot)
# saveWidget(PheWasPlotly, "PsyPheWasPlotly.html", selfcontained = T, libdir = "lib")